<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demje2</title>
  <meta name="keywords" content="demje2">
  <meta name="description" content="DEMJE2 Demonstrate nonlinear time series joint estimation for Mackey-Glass chaotic time series">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">examples</a> &gt; <a href="#">joint_estimation</a> &gt; demje2.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../menu.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\examples\joint_estimation&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>demje2
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>DEMJE2 Demonstrate nonlinear time series joint estimation for Mackey-Glass chaotic time series</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> DEMJE2 Demonstrate nonlinear time series joint estimation for Mackey-Glass chaotic time series

  The Mackey-Glass time-delay differential equation is defined by

            dx(t)/dt = 0.2x(t-tau)/(1+x(t-tau)^10) - 0.1x(t)

  When x(0) = 1.2 and tau = 17, we have a non-periodic and non-convergent time series that
  is very sensitive to initial conditions. (We assume x(t) = 0 when t &lt; 0.)

  We assume that the chaotic time series is generated with by a nonlinear autoregressive
  model where the nonlinear functional unit is a feedforward neural network. We use a
  tap length of 6 and a 6-4-1 MLP neural network (using the Netlab toolkit) with hyperbolic
  tangent activation functions in the hidden layer and a linear output activation.

   See also
   GSSM_MACKEY_GLASS
   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>	ADDRELPATH  Add a relative path which gets expanded into a absolute path</li><li><a href="../../.././ReBEL-0.2.7/core/datamat.html" class="code" title="function dm=datamat(x,M)">datamat</a>	DATAMAT  Packs a vector of data (length N) into a data matrix of dimension M-by-(N-M+1)</li><li><a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>	EKF  Extended Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>	GENINFDS  Generate inference data structure from a generalized state space model and user defined inference parameters.</li><li><a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>	GENSYSNOISEDS  Generate process and observation noise data structures for a given InferenceDS data structure</li><li><a href="../../.././ReBEL-0.2.7/core/mlpweightinit.html" class="code" title="function [W1, B1, W2, B2, W3, B3, W4, B4] = mlpweightinit(nodes)">mlpweightinit</a>	MLPWEIGHTINIT   Initializes the weights of a ReBEL MLP feedforward neural network</li><li><a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>	REMRELPATH  Remove a relative path (which gets expanded into a absolute path)</li><li><a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>	SRCDKF  Square Root Central Difference Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_mackey_glass.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_mackey_glass</a>	GSSM_MACKEY_GLASS  Generalized state space model for Mackey-Glass chaotic time series</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
</ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">% DEMJE2 Demonstrate nonlinear time series joint estimation for Mackey-Glass chaotic time series</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%  The Mackey-Glass time-delay differential equation is defined by</span>
0004 <span class="comment">%</span>
0005 <span class="comment">%            dx(t)/dt = 0.2x(t-tau)/(1+x(t-tau)^10) - 0.1x(t)</span>
0006 <span class="comment">%</span>
0007 <span class="comment">%  When x(0) = 1.2 and tau = 17, we have a non-periodic and non-convergent time series that</span>
0008 <span class="comment">%  is very sensitive to initial conditions. (We assume x(t) = 0 when t &lt; 0.)</span>
0009 <span class="comment">%</span>
0010 <span class="comment">%  We assume that the chaotic time series is generated with by a nonlinear autoregressive</span>
0011 <span class="comment">%  model where the nonlinear functional unit is a feedforward neural network. We use a</span>
0012 <span class="comment">%  tap length of 6 and a 6-4-1 MLP neural network (using the Netlab toolkit) with hyperbolic</span>
0013 <span class="comment">%  tangent activation functions in the hidden layer and a linear output activation.</span>
0014 <span class="comment">%</span>
0015 <span class="comment">%   See also</span>
0016 <span class="comment">%   GSSM_MACKEY_GLASS</span>
0017 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0018 <span class="comment">%</span>
0019 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0020 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0021 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0022 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0023 <span class="comment">%</span>
0024 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0025 <span class="comment">%   detail.</span>
0026 
0027 <span class="comment">%===============================================================================================</span>
0028 
0029 clc;
0030 clear all; close all;
0031 
0032 fprintf(<span class="string">'\nDEMJE2:  This demonstration shows how the ReBEL toolkit is used for joint estimation\n'</span>);
0033 fprintf(<span class="string">'         on a nonlinear time series (Mackey-Glass-30) problem. The scalar observation\n'</span>);
0034 fprintf(<span class="string">'         is corrupted by additive white Gaussian noise. A neural network is used as a\n'</span>);
0035 fprintf(<span class="string">'         generative model for the time series. We estimate both the model parameters and\n'</span>);
0036 fprintf(<span class="string">'         the underlying clean state from the noisy observations.\n'</span>);
0037 fprintf(<span class="string">'         We compare the performance of an EKF and a SRCDKF by iterating on the same sequence.\n\n'</span>);
0038 fprintf(<span class="string">'    NOTE : This demos is quite computationally expensive... so on a slow computer it might take a while.\n\n'</span>);
0039 
0040 
0041 <span class="comment">%--- General setup</span>
0042 
0043 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../gssm'</span>);         <span class="comment">% add relative search path to example GSSM files to MATLABPATH</span>
0044 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../data'</span>);         <span class="comment">% add relative search path to example data files to MATLABPATH</span>
0045 
0046 <span class="comment">%--- Initialise GSSM model from external system description script.</span>
0047 model = <a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_mackey_glass.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_mackey_glass</a>(<span class="string">'init'</span>);
0048 
0049 
0050 <span class="comment">%--- Load normalized Mackey glass data set</span>
0051 
0052 load(<span class="string">'mg30_normalized.mat'</span>);                           <span class="comment">% loads mg30_data from ../data/mg30_normalized.mat</span>
0053 
0054 mg30_data = mg30_data(100:100+300-1);
0055 
0056 
0057 <span class="comment">%--- Build state space data matrix of input data</span>
0058 
0059 X = <a href="../../.././ReBEL-0.2.7/core/datamat.html" class="code" title="function dm=datamat(x,M)">datamat</a>(mg30_data, model.statedim);                 <span class="comment">% pack vector of data into datamtrix for NN input</span>
0060 
0061 [dim,N]  = size(X);                                     <span class="comment">% dimension and number of datapoints</span>
0062 y  = zeros(model.obsdim,N);                             <span class="comment">% observation data buffer</span>
0063 
0064 clean_signal_var = var(mg30_data);                      <span class="comment">% determine variance of clean time series</span>
0065 
0066 SNR = 3;                                                <span class="comment">% 3db SNR</span>
0067 onoise_var = clean_signal_var/10^(SNR/10);              <span class="comment">% determine needed observation noise variance for a given SNR</span>
0068 
0069 model.oNoise.cov = onoise_var;                            <span class="comment">% set observation noise covariance</span>
0070 
0071 onoise = model.oNoise.sample( model.oNoise, N);   <span class="comment">% generate observation noise</span>
0072 
0073 y   = model.hfun( model, X, onoise);   <span class="comment">% generate observed time series (corrupted with observation noise)</span>
0074 
0075 <span class="comment">%----</span>
0076 
0077 ftype1 = <span class="string">'ekf'</span>;
0078 ftype2 = <span class="string">'srcdkf'</span>;
0079 
0080 
0081 <span class="comment">%--- Setup argument data structure which serves as input to</span>
0082 <span class="comment">%--- the 'geninfds' function. This function generates the InferenceDS and</span>
0083 <span class="comment">%--- SystemNoiseDS data structures which are needed by all inference algorithms</span>
0084 <span class="comment">%--- in the PiLab toolkit.</span>
0085 
0086 Arg.type = <span class="string">'joint'</span>;                                  <span class="comment">% inference type (state estimation)</span>
0087 Arg.tag = <span class="string">'Joint estimation for GSSM_MACKEY_GLASS system.'</span>;  <span class="comment">% arbitrary ID tag</span>
0088 Arg.model = model;                                   <span class="comment">% GSSM data structure of external system</span>
0089 
0090 InfDS = <a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>(Arg);                               <span class="comment">% create inference data structure</span>
0091 
0092 [pNoise1, oNoise1, InfDS1] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS,ftype1);    <span class="comment">% generate process and observation noise sources for EKF</span>
0093 [pNoise2, oNoise2, InfDS2] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS,ftype2);    <span class="comment">% generate process and observation noise sources for SRCDKF</span>
0094 
0095 
0096 <span class="comment">%--- Setup runtime buffers</span>
0097 
0098 Xh = zeros(InfDS.statedim,N);          <span class="comment">% state estimation buffer</span>
0099 Px = eye(InfDS.statedim);            <span class="comment">% initial state covariance</span>
0100 Px(model.statedim+1:<span class="keyword">end</span>,model.statedim+1:end) = 0.1*eye(model.paramdim);
0101 
0102 Xh(model.statedim+1:<span class="keyword">end</span>,1) = <a href="../../.././ReBEL-0.2.7/core/mlpweightinit.html" class="code" title="function [W1, B1, W2, B2, W3, B3, W4, B4] = mlpweightinit(nodes)">mlpweightinit</a>(model.nodes);              <span class="comment">% randomize initial model parameters</span>
0103 
0104 Xh1 = Xh;
0105 Px1 = Px;
0106 
0107 Xh2 = Xh;
0108 Sx2 = chol(Px)';                     <span class="comment">% SRCDKF is a square-root algorithm and hence it operates on the Cholesky factor</span>
0109                                      <span class="comment">% of the covariance matrix</span>
0110 number_of_runs = 10;                  <span class="comment">% we will iterate over the data 'number_of_runs' times</span>
0111 
0112 mse1 = zeros(1,number_of_runs);       <span class="comment">% buffers to store the MSE of each runs estimate</span>
0113 mse2 = mse1;
0114 
0115 mse1(1) = mean((y(1,:)-X(1,:)).^2)/var(y(1,:));    <span class="comment">% initial MSE of noisy signal</span>
0116 mse2(1) = mse1(1);
0117 
0118 <span class="comment">%--- Setup process noise data structures for joint estimation</span>
0119 
0120   pNoiseAdaptMethod = <span class="string">'anneal'</span>;                                <span class="comment">% setup process noise adaptation method (improves convergence)</span>
0121   pNoiseAdaptParams = [0.995 1e-7];                            <span class="comment">% annealing factor = 0.95     annealing floor variance = 1e-8</span>
0122 
0123   pNoiseCov0 = 1e-4*eye(model.paramdim);
0124 
0125   pNoise1.adaptMethod = pNoiseAdaptMethod;
0126   pNoise1.adaptParams = pNoiseAdaptParams;
0127 
0128   pNoise2.adaptMethod = pNoiseAdaptMethod;
0129   pNoise2.adaptParams = pNoiseAdaptParams;
0130 
0131   pNoise1.cov(2:<span class="keyword">end</span>,2:end) = pNoiseCov0;         <span class="comment">% set initial variance of process noise parameter estimation subvector</span>
0132   pNoise2.cov(2:<span class="keyword">end</span>,2:end) = chol(pNoiseCov0)';  <span class="comment">% set initial variance of process noise parameter estimation subvector</span>
0133 
0134 
0135 <span class="comment">%---</span>
0136 
0137 fprintf(<span class="string">'\n Running joint estimators ... \n\n'</span>);
0138 
0139 
0140 <span class="comment">%--- Call inference algorithm / estimator</span>
0141 
0142 <span class="keyword">for</span> k=1:number_of_runs,
0143 
0144   fprintf(<span class="string">' [%d:%d] '</span>,k,number_of_runs);
0145 
0146 
0147   <span class="comment">%------------------- Extended Kalman Filter ------------------------------------</span>
0148 
0149 
0150   [Xh1, Px1, pNoise1] = <a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>(Xh1(:,1), Px1, pNoise1, oNoise1, y, [], [], InfDS1);
0151 
0152 
0153   <span class="comment">%------------------- Square-root Central Difference Kalman Filter -------------</span>
0154 
0155   InfDS2.spkfParams = sqrt(3); ;                                 <span class="comment">% scale factor (CDKF parameter)</span>
0156 
0157   [Xh2, Sx2, pNoise2] = <a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(Xh2(:,1), Sx2, pNoise2, oNoise2, y, [], [], InfDS2);
0158 
0159   <span class="comment">%---------------------------------------------------------------------------------</span>
0160 
0161 
0162   <span class="comment">%--- Calculate normalized mean square estimation error</span>
0163 
0164   mse1(k+1) = mean((Xh1(1,:)-X(1,:)).^2)/var(y(1,:));
0165   mse2(k+1) = mean((Xh2(1,:)-X(1,:)).^2)/var(y(1,:));
0166 
0167   <span class="comment">%--- Plot results</span>
0168 
0169   figure(1); clf; subplot(<span class="string">'position'</span>,[0.025 0.1 0.95 0.8]);
0170   p1 = plot(X(1,:),<span class="string">'b'</span>,<span class="string">'linewidth'</span>,2); hold on
0171   p2 = plot(y,<span class="string">'g+'</span>);
0172   p3 = plot(Xh1(1,:),<span class="string">'m'</span>);
0173   p4 = plot(Xh2(1,:),<span class="string">'r'</span>); hold off
0174   legend([p1 p2 p3 p4],<span class="string">'clean'</span>,<span class="string">'noisy'</span>,<span class="string">'EKF estimate'</span>,<span class="string">'SRCDKF estimate'</span>);
0175   xlabel(<span class="string">'time'</span>);
0176   ylabel(<span class="string">'x'</span>);
0177   title(<span class="string">'DEMSE3 : Mackey-Glass-30 Chaotic Time Series Joint Estimation'</span>);
0178 
0179   figure(2);
0180   p1 = plot(mse1(2:k+1),<span class="string">'m-o'</span>); hold on;
0181   p2 = plot(mse2(2:k+1),<span class="string">'r-s'</span>); hold off;
0182   legend([p1 p2],<span class="string">'EKF'</span>,<span class="string">'SRCDKF'</span>);
0183   title(<span class="string">'Normalized MSE of Estimates'</span>);
0184   xlabel(<span class="string">'k'</span>);
0185   ylabel(<span class="string">'MSE'</span>);
0186   drawnow
0187 
0188   fprintf(<span class="string">'  Mean-square-error (MSE) of estimates : EKF = %4.3f    SRCDKF = %4.3f\n'</span>, mse1(k+1), mse2(k+1));
0189 
0190 
0191   <span class="comment">%-- Copy last estimate of model parameters to initial buffer position for next iteration...</span>
0192 
0193   Xh1(model.statedim+1:<span class="keyword">end</span>,1) = Xh1(model.statedim+1:<span class="keyword">end</span>,end);              <span class="comment">% copy model parameters over</span>
0194   Xh1(1:model.statedim,1) = zeros(model.statedim,1);                        <span class="comment">% reset state estimate</span>
0195   Px1_temp = eye(InfDS.statedim);                                           <span class="comment">% copy covariance of parameter estimates</span>
0196   Px1_temp(model.statedim+1:<span class="keyword">end</span>,model.statedim+1:end) = Px1(model.statedim+1:<span class="keyword">end</span>,model.statedim+1:end);
0197   Px1 = Px1_temp;
0198 
0199   Xh2(model.statedim+1:<span class="keyword">end</span>,1) = Xh2(model.statedim+1:<span class="keyword">end</span>,end);              <span class="comment">% copy model parameters over</span>
0200   Xh2(1:model.statedim,1) = zeros(model.statedim,1);                        <span class="comment">% reset state estimate</span>
0201   Sx2_temp = eye(InfDS.statedim);                                           <span class="comment">% copy covariance of parameter estimates</span>
0202   Sx2_temp(model.statedim+1:<span class="keyword">end</span>,model.statedim+1:end) = Sx2(model.statedim+1:<span class="keyword">end</span>,model.statedim+1:end);
0203   Sx2 = Sx2_temp;
0204 
0205 
0206 <span class="keyword">end</span>
0207 
0208 
0209 <span class="comment">%--- House keeping</span>
0210 
0211 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../gssm'</span>);       <span class="comment">% remove relative search path to example GSSM files from MATLABPATH</span>
0212 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../data'</span>);       <span class="comment">% remove relative search path to example data files from MATLABPATH</span></pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>